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This paper compares a fluid/thermal simulation, in Fluent, with a low-g, nitrogen slosh 
and boiling experiment. In 2010, the French Space Agency, CNES, performed cryogenic 
nitrogen experiments in a low-g aircraft campaign. From one parabolic flight, a low-g 
interval was simulated that focuses on low-g motion of nitrogen liquid and vapor with 
significant condensation, evaporation, and boiling. The computational results are compared 
with high-speed video, pressure data, heat transfer, and temperature data from sensors on 
the axis of the cylindrically shaped tank. These experimental and computational results 
compare favorably. The initial temperature stratification is in good agreement, and the 
two-phase fluid motion is qualitatively captured. Temperature data is matched except that 
the temperature sensors are unable to capture fast temperature transients when the sensors 
move from wet to dry (liquid to vapor) operation. Pressure evolution is approximately 
captured, but condensation and evaporation rate modeling and prediction need further 
theoretical analysis. 
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Nomenclature 

acceleration, acceleration of center of gravity, m/s 2 
Air Liquide Advanced Technology 
Centre National d’Etudes Spatiales 

specific heat at constant pressure and constant volume, J/kg-K 

VOF fraction, unitless 

gravitational acceleration, m/s 2 

thermal conductivity, W/m-K 

molecular weight, kg/kmol 

mass flux due to evaporation (>0) and condensation, kg/s-m 2 
static pressure, Pa, saturation pressure at temperature, Pa 
universal gas constant, J/kmol-K 
radius vector from center of gravity, m 
temperature, K, saturation temperature at pressure, K 
User Defined Function (Fluent) 

Volume of Fluid 

fluid velocity, m/s, radial velocity relative to center of gravity, m/s 

dissipation per unit turbulent kinetic energy, 1/s 

compressibility factor, unitless 

time step, s 

surface tension, N/m 

turbulent kinetic energy, m 2 /s 2 

molecular viscosity, Pa s 

density, kg/m 3 
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cr C ond> <^eva P - accommodation coefficients for condensation and evaporation, unitless 
r = time constant for temperature sensor, s 

CO, CO = angular velocity vector, rad/s, angular acceleration vector, rad/s 2 

iiq, vap = subscripts indicating liquid and vapor phases 


I. Introduction 

I N 2010, Centre National d’ Etudes Spatiales and Air Liquide Advanced Technologies conducted a parabolic flight 
campaign which included an experimental test bench, CryOgenic, to study nitrogen sloshing and boiling in low-g 
conditions. Computational fluid dynamics simulations accompanied these experiments [1] [2]. This experimental 
and computational work is part of a long term effort by CNES to examine cryogenic fluid behavior in low-g 
conditions, including evaporation, condensation, boiling, stratification, pressurization and de-pressurization, as well 
as helium pressurization. 

CNES is interested [3] in cryogenic fluid behavior in flight conditions, gathering data for validating 
computational simulations, and predicting tank pressure, thermal stratification, and tank outlet/engine inlet 
conditions for engine restart after an orbital coast phase. NASA is interest in long term storage of cryogenic 
propellants for on-orbit refueling and long duration space missions. Cryogenic propellants promise higher specific 
impulse than storable hypergolic fuels, but storage for long duration missions must be demonstrated. A pacing 
technology for manned Mars missions using Nuclear Thermal Propulsion is storage of many tons of liquid hydrogen 
propellant for several years with minimal boil-off loss. With this convergence of interests, NASA and CNES agreed 
to benchmark simulations, and this is one of the test cases. 

Geometry and Grid 

The experimental apparatus consisted of nitrogen, gas and liquid, in a 
sapphire cylindrical shell with aluminum and stainless steel (inox) 
lids (Figure 1) all contained in an insulating vacuum chamber. To 
remove the small residual heat fluxes into the vessel (and reduce boil 
off), the lower lid contains a cryocooler. Instrumentation included 
high speed video and twelve temperature sensors mounted to a post 
along the axis of the shell. The dimensions of the cylindrical 
nitrogen vessel are approximately 6 cm by 10 cm, with a slosh 
frequency of near 4 Hz. 

The three-dimensional grid is shown in Figure 2. It is a full 360- 
degree sector of the apparatus with 569,110 fluid cells, and 685,858 
solid cells. The solid grid is unstructured with variable resolution. 
The interior of the fluid grid is a uniform, structured grid with -1 mm 
resolution, but near the fluid- solid interface, it becomes unstructured 
and merges into the solid grid. The grid does not include the central 
post where the fluid temperature sensors are mounted. The region of 
the joints and sealing gaskets between the sapphire shell and metal 
lids is meshed; however, thermal isolation is enforced along the 
interface lines shown in Figure 2. The grid was partitioned to run of 
16 or 32 processors. 

III. Numerical Methods and Fluent Settings 

ANSYS Fluent version 13 [4] was used to solve thermal equations in the solid coupled to thermal/fluid equations 
in the fluid region. While the energy equation is solved in the solid region, mass, momentum, energy, and 
turbulence equations are solved in the fluid region with second-order upwind schemes. The PISO scheme was used 
for the pressure-velocity coupling, and the PRESTO ! scheme was used for the pressure interpolation. The gas phase 
was modeled as an ideal gas, and a Boussinesq approximation was used for the liquid phase. Two-phase flow was 
resolved using the Volume of Fluid method [5]. The k-w SST turbulence model of Menter [6] [7] was used with a 
turbulent damping coefficient of 10. The simulation was time-accurate (second order implicit temporal scheme) 
with a time step At = l.OxlO" 4 seconds. The simulation was run on Pleiades at the NASA Advanced 
Supercomputing Facility. Add execution time. 
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II. 



Figure 1. Cross-sectional geometry for 
experimental apparatus. Picture from 
Mathey et at. [1]. 



Figure 2. Grid cross-section for fluid and solid regions of experimental apparatus. Grid is 3 -dimensional and a full 360 

degree sector. 


The following sections explain the fluid and material properties used, mass and heat transfer due to evaporation 
and condensation, aircraft low-g acceleration, and thermal, boundary, and initial conditions. 

A. Temperature and Pressure Dependence of Thermophysical Properties 

Constant fluid physical properties at reference conditions (77.224 K, 10 5 Pa) are in error [8] by as much as 10% 
at extreme simulation conditions (110 K, 3x1 0 5 Pa). However, linear temperature variations accurately represent 
viscosity, //, thermal conductivity, k , surface tension, y, and the liquid’s specific heat, C p . There are both 
temperature and pressure variations for the gas’s specific heat, C p , the heats of vaporization and condensation. The 
compressibility factor, Z, varies about 8% over the temperatures and pressures in the experiment. Z was not 
corrected (ideal gas assumption), but temperature variations were included for the gas’s C p and the heats of 
vaporization and condensation. 

For the solid materials, temperature dependent data (specific heat, C p , and thermal conductivity, k ) were provided 
by CNES [9], and 310 stainless steel agrees with this data. But NIST [10] and MMPDS [11] data do not agree with 
the aluminum thermal conductivity. 

The saturation line conditions for nitrogen, P sat (T) and T sat (P), were fit to Reynolds [12] and NIST [8] data, 
respectively. This is a significant improvement over the Clausius -Clapeyron equation, particularly away from 
reference conditions. 


B. Mass Transfer and Heat of Vaporization / Condensation 

The rate of mass transfer and heats of vaporization/condensation are calculated from the Schrage equation, Eq. 
(1), which is derived, in turn, from Maxwell’s distribution. This evaporation/condensation internal boundary 
condition enforces saturation conditions at the liquid/vapor interface, and it is implemented in a Fluent User- 
Defined-Function (UDF). 

Equation (1) is modified by assuming a cond = cr evap ; T vap = T Uq ; P vap = P sat . Implementation requires a length 
scale, sqrt(l/lgrad cl 2 ), where c is VOF fraction, to convert the mass flux (kg/s-m 2 ) to a volumetric rate, that is, to 
distribute the mass flux over the width of the liquid/vapor interface [13]. 
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Unfortunately, Eq. (1) and its implementation are not predictive of mass transfer due to evaporation or 
condensation across an interface in this simulation, and results are matched by adjusting the accommodation 
coefficient, cr. The best-fit values of the accommodation coefficient are cr = 1.0x1 O' 4 for evaporation and 
condensation at the liquid/vapor interface and for vapor phase condensation (small in this simulation). For boiling 
(evaporation initiated away from the liquid/vapor interface), which dominates in this simulation, cr = 5.0xl0‘ 3 was 
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used. This boiling model includes a superheat criterion, T max - T sat (P) > 5 K, that is, the liquid temperature within 
the cell must exceed the saturation temperature by a fixed temperature in the liquid phase. T max is the maximum 
temperature in the cell, including adjacent walls. The vapor phase condensation model includes a similar 
supercooled liquid criteria based on the minimum temperature. 

Data exists for boiling heat transfer rates for 1 -g, nitrogen boiling in wires [14] and cylinders [15]; however, this 
wire and cylinder data is not clearly applicable here, so they can only provide guidance. 


C. Non-Inertial Reference Frame Accounts for Low-g Aircraft Acceleration 

A non-inertial reference frame will 
account for the acceleration of the aircraft 
during low-g (and high-g) parabolic flight. 

In the general case, such as a satellite with 

measured linear acceleration, a cg , angular 


velocity, (D , and angular acceleration, (b , 
the acceleration components are given by 
Eq. (2). The terms included on the RHS of 
the x-, y-, and z-momentum equations are 

p a , and in the energy equation, pa V . 

In this experiment, angular velocity and 
acceleration are zero. Only two 
components of the linear acceleration of the 
aircraft, a x and a z , were measured, and a y is 
assumed to be zero. The acceleration 
components were measured during 
parabolic flight at a 10 Hz sampling rate, as 
shown in Figure 3. Values are interpolated 
as piece- wise linear fits. The steady initial 



Figure 3. Acceleration profile during low-g interval. High-g refers 
to the preceding and following intervals. 


conditions assumed in the high-g phase of parabolic flight are a x = -16.5 m/s and a z = - 1.93 m/s . 


a=5 +2^xv +fflxr + fi}x(fflx r) 

eg r 


( 2 ) 


The Bond or Eotvos number is in the range [0.3, 6.] 
during the low-g phase, and it is approximately 1 for 
about 10 seconds. 

The weightless interval of parabolic flight is referred 
to as low-g, while high-g refers to the intervening 
interval. 

D. Thermal, Boundary, and Initial Conditions 

Despite being contained in a vacuum chamber, small 
heat fluxes exist into the experimental apparatus, and 
they are compensated by a cryocooler attached to the 
lower aluminum lid. The surface heat fluxes are due to 
radiation from the vacuum chamber inner surface, and 
conduction through the near vacuum conditions. 
Consequently, heat flows through the apparatus as 
discussed in the next section. The cryocooler removes 
about 4 W, and the surface heat inflow is estimated and 
assumed to be uniformly distributed over the top lid (1 
W), bottom lid (1 W), and the sapphire surface (2 W). 

Nitrogen fills the tank 55% full. The liquid-to- 
vapor contact angle is taken as 5 degrees. Between low- 
g intervals in the parabolic flight, the conditions are 
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Figure 4. Initial temperature determined from a 
transient simulation of the high-g interval. Temperature 
sensor locations are also marked. Liquid/vapor interface 
indicated between sensors tl2c and tl2e. High temperature 
gradient near tl2b. 
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assumed to be steady, yet since a z * 0, the liquid/vapor interface is tilted. This interface position and the initial 
thermal conditions are calculated with a transient analysis, and a steady thermal/fluid solution was achieved well 
within the 90 seconds of physical time simulated. This initial temperature distribution is shown in Figure 4. This 
initial solution was extended into the transient, low-g simulation, that is, the non-inertial reference frame UDF 
commenced with non-constant values of a x and a z . 

TV. Operational Behavior of the Experimental Apparatus 
E. How Does Boiling Occur in Low-g Conditions? 

Thermal isolation of the sapphire shell from the metal lids requires a large temperature difference for heat to 
diffuse through the relatively non-conductive nitrogen gas during the high-g interval. Consequently, a hot, top lid 
develops during the high-g interval (Figure 4), and boiling occurs during low-g conditions when liquid nitrogen re- 
orients and impinges on this hot lid. 

A thermal analysis of the apparatus confirms this scenario for boiling in low-g conditions. Small heat fluxes 
enter the external surfaces of the sapphire shell and lids, and this heat is removed by the cryocooler. What path will 
this heat take through the solid walls or the nitrogen gas and liquid? The thermal conductivity of the sapphire shell 
and aluminum and stainless steel lids is at least 3 orders of magnitude greater than the nitrogen vapor’s conductivity 
(Figure 5-left), and two orders of magnitude greater than nitrogen liquid [9] [8] [12]. A spreadsheet analysis shows 
that if either lid thermally contacts the sapphire shell, the heat should quickly diffuse through the solid walls with a 
temperature difference of less than a degree. 



Figure 5. Relative thermal conductivity of materials (left), and flow of heat during high-g conditions (right). The 

sapphire cylinder and metal lids have thermal conductivities 3 orders of magnitude above the nitrogen gas; with thermal 
isolation at the joints , heat must diffuse through the gas, a large temperature difference develops, and a hot top lid. 

Seals 3 * 5 between the sapphire shell and metal lids prevent nitrogen loss during operation, but they are also 
thermally insulating. Hence, heat is forced through the nitrogen vapor and liquid, and a much higher temperature 
gradient is required, and the top, stainless lid can be 30 K hotter than the sapphire and aluminum — enough to boil 
liquid nitrogen when it re-orients during low-g conditions. 

Consequently, thermal isolation of the sapphire shell and the metal lids is essential to simulating this experiment. 

F. What is the Experiment’s Fluid/Thermal Behavior? 

The experiment’s behavior is complex. During high-g conditions, liquid nitrogen settles to the bottom (Figure 
4), the fluid is relatively motionless, and a hot lid develops as explained in Section E. During low-g conditions, 
acceleration reverses sign several times (Figure 3) and the liquid re-orients (i.e. from bottom-to-top) with each 

3 With large temperature changes and different thermal expansion coefficients, direct contact between sapphire and 
metal lids is hard to control, and the default design is a flexible seal between the mating surfaces and use of 
Belleville washers which, like springs, maintain a pre-load tension while accommodating thermal deflections. 
Bintley et al. [16] measured thermal conductance through sapphire- sapphire and sapphire-metal bolted joints and 

found good thermal isolation. 
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acceleration sign change. Since the lid temperature is ~30 K above the fluid saturation temperature, T sat (P), rapid 
boiling occurs when liquid nitrogen contacts the lid. Bubbles do not move away from the lid (due to buoyancy) as 
quickly as in 1 -g conditions, so larger bubbles are expected. As the hot gas bubbles move through the relatively 
cold liquid, there is significant vapor condensation near the liquid/vapor interface, bubbles decrease in size 
noticeably, and heat of condensation is transferred to the liquid. This condensation tends to balance boiling, reduces 
the pressure increase, and quickly moves large amounts of heat through the liquid/vapor and the apparatus. 

V. Results and Comparison with Experimental Data 

The computational results were compared with the experimental data in five ways: initial thermal profile, high- 
speed video, pressure evolution data, net heat transfer, and temperature probe data. 

G. Comparison with Initial Thermal Profile 

As noted in Section E, the boiling Table 1. Comparison of experimental data and the initial thermal profile 

hot lid is established during high-g from a steady-state simulation of high-g conditions. Temperature sensor 

intervals in the parabolic flight as the locations are shown in Figure 4. Colors indicate solid, liquid, and gas. 
small heat flux is forced through the 
low thermal conductivity gas phase. 

This steady-state condition was 
simulated in a transient calculation 
and used as a starting point for 
simulation of the low-g interval. 

Fluid temperature sensor data is 
available for comparison, but wall 
temperature data is not. The 
comparison is shown in Table 1. 

The agreement is very good 
between experimental and 
computational results. The 

disagreement is largest in the high 
temperature gradient (tl2b) near 
where the top lid meets the sapphire 
window at a seal. A smaller 
discrepancy exists near the bottom 
wall, and it may be due to the 
estimated thermal boundary 
conditions (Section D) or the 
aluminum thermal conductivity. 

H. Visual Comparison with High-Speed Video 

The high-speed video compares well with the simulation results. During the initial re-orientation, the 
deformation of the liquid/vapor interface (and timing) is captured well, as shown in Figure 6. The liquid/solid 
contact line is also captured as shown by the arrows in Figure 6. 

The liquid/vapor experiences a number of re -orientations as the vertical acceleration changes sign (Figure 3), and 
the timing of these events is well captured. The speed of the bubbles is visually similar as acceleration changes. 
The simulation also captures the initiation of boiling on the top lid. As would be expected, the bubbles ’ size and 
shape are not duplicated exactly, but they are captured qualitatively (Figure 7). However, the simulation appears to 
have fewer bubbles than the experimental movie. The bubbles in the movie tend to be more elongated (higher 
aspect ratio) than those in the experiment. Experimentally, the bubbles have large, chaotic surface waves, while 
computationally they are relatively smooth; Figure 7 illustrates this difference. The simulation does capture the 
merging of bubbles. 


Sensor Location 

CNES 

Measured 

Temperature 

(K) 

Steady 

Simulation 

Temperature 

(K) 

Percent 

Difference 

(%) 

Top Lid Center 


105.1 


Top Lid Edge 


104.0 


Top Lid Side 


103.3 


tl2a 

103.77 

104.1 

0.3 

tl2b 

86.00 

88.9 

3.4 

tl2c 

77.41 

78.1 

0.9 

tl2d 

77.24 

77.9 

0.9 

tl2e 

77.21 

77.2 

0.0 

tl2f 

77.17 

77.1 

-0.1 

tl2g 

77.14 

77.0 

-0.2 

tl2h 

77.09 

77.0 

-0.2 

tl2i 

77.02 

76.9 

-0.2 

tl2j 

76.93 

76.7 

-0.3 

tl2k 

76.42 

76.3 

-0.1 

tl21 

76.00 

76.3 

0.4 

Bottom Lid Center 


75.3 
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Figure 6. Comparison of computational predictions (left) and experimental data (right) showing the liquid/vapor 
interface as it rapidly deforms with the initial reduction in gravity. Arrows indicate the liquid/solid contact line. At left, the 

blue bands indicate the sapphire/lid joints. 



Figure 7. Comparison of computational predictions (left) and experimental data (right) showing the liquid/vapor interface 
during the intense boiling phase. At left, the blue bands indicate the sapphire Aid joints. 


I. Pressure Data 

Pressure is a measure of the net mass transfer due to evaporation and condensation. Figure 8 (left) shows the 
time evolution of mass transfer rate (evaporation, condensation, and net), and (right) shows pressure with time for 
two accommodation coefficients. Evaporation (boiling) precedes condensation as the liquid nitrogen contacts the 
top lid. After a surge of condensation and a pressure drop, evaporation and condensation reach an approximate 
balance with no apparent correlation with subsequent re -orientations. 

The computational results for pressure, shown in Figure 8, show general agreement with the experimentally 
measured data. The simulation does capture the initial rapid pressure rise and reversal, but it does not capture the 
gentle pressure rise between 97 and 107 seconds. 

J. Heat Transfer Rate and Fluid Enthalpy Increase 

We verify the boiling heat transfer rate by calculating the increase in fluid enthalpy in the simulation and 
comparing with an estimate from the experiment. The results, shown in Figure 9, show good agreement, and 
indicate the accommodation coefficient value, cj = 5.0xl0‘ 3 , used in Eq. (1) is a good fit. The computational 
estimate evaluated the integral of Eq. (3), and the experimental estimate of enthalpy assumed the fluid was thermally 
stratified at the temperatures indicated by the probes before and after the low-g interval. 
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Figure 8. Evaporation, condensation, and net mass transfer rates (left), and pressure evolution for two accommodation 
coefficients compared with experimental measurements (right). The pale red line indicates acceleration, a/g. 


f V0 lPliq(T 1 ref ) ^v iiq VO^liq dVol 


( 3 ) 


K. Comparison of Temperature Probe Results 

Although experimental wall temperature data were not available, predictions from the simulation are given in 
Figure 10. From 114 to 120 seconds in Figure 10, the top lid is reheating. How long does the lid take to reheat 
during the high-g interval? Combining Figure 10 data with the transient solution for initial conditions, the lid re- 
heating time is expected to be completely reheated within 60 seconds. 

The temperature evolution at two temperature probes is shown in Figure 11, and these are typical of the 12 
probes. The results for probe tl2a show the probe is in the hottest vapor at the top of the tank. As the hot vapor 
moves away and is replaced with cold liquid, the predicted quick temperature drop agrees with the experimental 
measurements. From 95 to 113 seconds, the liquid slowly heats, and there is general agreement between the two 
results. However after 115 seconds, the computation predicts re -heating of the vapor near the lid, while the 
experimental data do not show this. One explanation is that liquid nitrogen remains in the fill lines above probe 
tl2a, and this liquid drains over the temperature probe after the re -orientation. The fill lines are not modeled in the 



Time (s) 

Figure 9. To check if the heat transfer rate is correct, compare the computationally measured change in fluid enthalpy 
with an estimate from the experiment. The pale red line indicates acceleration, a/g. 
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Figure 10. Time evolution of predicted wall temperature for four probes in the computational simulation. Probe 
locations are indicated by stars in the graphic at right. The pale red line indicates acceleration, a/g. 

During the initial re -orientation of the fluid (~93 s), the hot gas at the lid moves to the bottom of the tank, passing 
the temperature probes mounted on the axis, and liquid nitrogen occupies the top of the tank, contacts the lid, and 
begins to boil. The movement of this hot gas is clear in the computational video, and in the computational 
temperature predictions for probe tl2g, shown in Figure 11 (right) between 94 and 95 seconds. However, it is not 
measured by probe tl2g. In fact, while probes tl2a and tl2b capture the dry-to-wet temperature jump, probes tl2d 
to tl21 do not capture the wet-to-dry temperature jump as the hot gas moves past the temperature sensors. 

Typical diode temperature sensors have a time constant, r, near 0.1s [16]. Consequently 95% of a temperature 
jump is captured in three time constants, or 0.3 s. The computationally predicted hot vapor exposure time of the 
temperature probes is 0.3 to 0.5s, hence the sensors could detect most of the jump. However, delays in wet-to-dry 
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Figure 11. Predicted and measured temperatures at sensors tl2a (left) and tl2g (right). The pale red line indicates 

acceleration, a/g. 
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response time have been observed experimentally and analyzed [17]. The explanation is that wet-to-dry transition of 
a diode sensor includes a liquid film that must vaporize before gas temperature is measured. 


VI. Conclusion 

The computational and experimentally results agree. The initial temperature profile, bubble motion, and heat 
transfer, are in good agreement, while the pressure evolution, and shape and number of bubbles is only acceptable. 
Evaporation/condensation model is not predictive of rates and more theoretical analysis is needed. Turbulence 
damping keeps turbulence away from interface. 
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